
obs_success <- 12
n <-20

alpha_prior <- 3.06
beta_prior <- 2.56

poster_alpha <- alpha_prior + obs_success
poster_beta <- beta_prior + n - obs_success

n_sim <- 20000

theta_prior <- rbeta(n_sim, alpha_prior, beta_prior)
y_prior <- rbinom(n_sim, n, theta_prior)

summary(y_prior)

prob_prior <- mean(y_prior >= obs_success)
cat("Prior probability of observing at least", obs_success, "successes:", prob_prior, "\n")

hist(y_prior, breaks = 20, main = "Prior Distribution of Successes", 
     xlab = "Number of Successes", col = "lightblue", probability = TRUE)

abline(v = obs_success, col = "red", lwd = 2)


theta_post <- rbeta(n_sim, poster_alpha, poster_beta)
y_post <- rbinom(n_sim, n, theta_post)

prob_post <- mean(y_post >= obs_success)
cat("Posterior probability of observing at least", obs_success, "successes:", prob_post, "\n")

hist(y_post, breaks = 30, main = "Posterior Distribution of Successes", 
     xlab = "Number of Successes", col = "lightgreen",
     probaility = TRUE)

abline(v = obs_success, col = "red", lwd = 2)

curve(dbeta(x, alpha_prior, beta_prior), 
      from = 0, to = 1, add = TRUE, col = "red", lwd = 2)
curve(dbeta(x, poster_alpha, poster_beta), 
      from = 0, to = 1, add = TRUE, col = "blue", lwd = 2)

# legend("topright", legend = c("Prior", "Posterior"), 
#        col = c("red", "blue"), lwd = 2)


beta_intergrate <- function(a, b) {
  integrate(function(x) dbeta(x, a, b), 0, 1)$value
}

binomial_integrate <-function(n, k, a, b) {
  integrate(function(x) dbinom(k, n, x) * dbeta(x, a, b), 0, 1)$value
}

m <- binomial_integrate(n, 1, alpha_prior, beta_prior)

# integrate for x in 1 and n
results <- sapply(1:n, function(k) binomial_integrate(n, k, alpha_prior, beta_prior))
results <- c(0, results)

summary(results)
sum(results)
mean(results)
sd(results)

lines(results, type='l', col='red', lwd=2)

plot(results, type='l', color='red', lwd=2)


barplot(results, names.arg = 0:n, 
        main = "Expected Number of Successes with Prior", 
        xlab = "Number of Successes", ylab = "Probability",
        col = "lightblue")


xm <- choose(35,31)
m <- integrate(function(x) xm * x^31*(1-x)^4, 0, 1)$value
